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1. Introduction 

A crucial problem of event generation at next-to-leading order (NLO) in perturbation theory 
originates from the loop contributions to multi-parton matrix elements. The number of individual 
Feynman graphs rises dramatically with the number of external legs, and tensor reduction methods 
increase the number of terms even more. This poor scaling behavior, and the associated instabilities 
due to large numerical cancellations in matrix elements, make it worthwhile to investigate different 
and alternative approaches to the problem of automatic NLO computation and simulation. 

We propose a new method for the evaluation of loop integrals, which allows for direct nu- 
merical computation without a decomposition into a set of basic integrals and corresponding large 
coefficients. 1 These proceedings summarize the work of a recent publication [Qj. Here, the matrix 
element is re-expressed using an improved version of the Feynman Tree Theorem (FTT) [Qj. The 
loop integrals are transformed into phase space integrals and are evaluated numerically along with 
the phase space integral over the external partons of the process under consideration. Thus, there 
is only one step of numerical integration involved in the computation of any particular integrated 
cross section or distribution. We therefore can make use of powerful existing technologies for 
numerical phase-space integration by tree-level event generators to evaluate NLO processes. 



2. Feynman Tree Theorem 

The integrand I(k) of a loop integral over momentum k can be written as a product of Feynman 
propagators F times a regular function N(k) in the numerator. Suppressing possible indices in the 
numerator, we have 

I{k)=N(k)]jF(k + p i ,m i ). (2.1) 

The pi are linear combinations of external momenta, m, the masses of the physical particle the 
propagator corresponds to. We define 

Fi = F(k + pi,mi) = — — Y—, (2.2) 

(k + pi) 1 - mf + is 

using 't Hooft-Feynman gauge throughout this discussion. 

The Feynman Tree Theorem [||] states, that the loop integral over (2.1) can be replaced by 
an integral over a sum of terms, where in each term one or more propagators F, are replaced by 
8-f unctions, 



J dk-i{k) =j dk-N(k) a'f • • • -£ a'a'f ••• + ...- (-i)"£a' • • • a' 



(2.3) 



with dk = d 4 k • (2^)" 4 . The 5-functions A\ = j§-8(k° - (-/>?+£,-)) cancel the ^-integration and 
cut the loop by setting the momentum of the original propagator Fj on its mass-shell. The numer- 
ator structure of this propagator factorizes into a product of wave functions of the corresponding 
particle. The original loop integral can therefore be written as a sum of tree graphs, with a phase 
space integration replacing the loop integral. 



'This is in contrast to methods used for recent computations of 2 —¥ 4 processes at the LHC These either use 
Passarino-Veltman style reductions or a unitarity based approach to decompose amplitudes into a set of basic integrals. 
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Derivations of the FTT in its form (2.3) can also be found in [Q] and [J3||. In the latter work, 
Catani et al. exploited the relation between loop integrals and phase space integrals to express loop 
amplitudes as sum of tree amplitudes arising solely from single cuts. Contributions from multiple 
cuts, which are present in (|2j|), are compensated by a non-trivial ie -prescription in the propagators. 

To make use of the FTT in a numerical evaluation of loop integrals involving only real num- 
bers, at some point we have to set the ie terms in the denominators to zero. Using the identity 

1 @— — =F«r5(jt-a), (2.4) 



x — a ± ie 



where & is Cauchy's Principal Value, we can rewrite (2.2) after a partial fraction decomposition: 



F i = P i +-A l i + ~A u i , (2.5) 

where we defined A" = jE.S(k° — (—p® —£";)), as the delta function that sets the associated four- 
momentum on its mass-shell with negative zero component. Pj stands for the propagator with no 
ie -prescription in the denominator. 



Replacing the propagators F; in (2.3) by (2.5), we obtain a version of the FTT, which is better 



+ 



suited for a direct numerical integration: 

J dk ■ I(k) =Jdk-N(k) [A[P 2 ■■■P n +P l A l 2 P 3 ---P n + ...+P r -- Pn-iK] 

~jdk-N(k) £ JL(l-(-l)*) AW. (2.6) 

perm. 
U + L>2 

The sum runs over all possible permutations, where the functions (A l ,A",P) appear (L,U,P) times. 

We explicitly wrote out the terms containing one A' function in the first line of ( |2.6| ). After 
k° -integration, all of these terms can be interpreted as tree graphs with one additional incoming 
and outgoing particle and an additional phase space integral over this particles momentum. This 
integral can now be pulled out of the individual graphs and put in front of the amplitude, alongside 
the phase space integrals over the external particles which are present in a cross-section calculation. 
Since the extra integration is also of the form of a phase-space integral, techniques developed for 
the integration, in particular multi-channel sampling, can immediately be adopted. Furthermore, 
this integration can be performed simultaneously with the external phase space integration, which 
is in contrast to common analytical methods, where for each individual configuration of external 
momenta the analytical result of the loop integration has to be numerically evaluated. 

If the momenta of two or more propagators go on-shell simultaneously in the momentum inte- 



gration, terms with more than one A-function in (2.6) are non- vanishing. These give rise to further 



contribution to the final result, but can be evaluated rather easily, since after the &° -integration, the 
remaining 8-f unctions lower the dimension of integration to two or less. 

3. UV and IR Subtraction 

As long as we are dealing with massive theories, a convenient set of renormalization condi- 
tions is given by the on-shell renormalization scheme. This requires the pole of the real part of 
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propagators to coincide with the corresponding particle mass, with residue 1. Vertex functions are 
set equal to the tree level vertex functions for on-shell external legs. 

We make use of a variation of the BPHZ procedure, which results in loop graphs acting as 
counterterms. These can be evaluated under the same phase space integral over the loop momen- 
tum. UV divergences are then cancelled locally in three-momentum space. 

Consider a 1PI one-loop graph T n (p\, . . . ,p n ) with superficial degree of divergence coiT). We 
define the T operator as a Taylor expansion around on-shell momenta pj, with pj = mj: 



n-l npn 

t o r>! ,..., Pn ) = , . . . ,p n ) + £ ( Pi - piY 

i 



+ 

Pl=Pl,-,Pn=Pn 



dp* 

.+ (3.1) 



i 7i- 1 ^dyn 

h,-Jd 



(l - i i, d Pu--- d Pi 



p l= p u ...,p„=p n 



up to d = co(r). With this T operator, the renormalized 1PI n-point functions 

r n ( Pl , Pn ) = r n ( Pl , Pn ) - ToT n ( Pu . . . , Pn ) (3.2) 

fulfill the renormalization conditions of the on-shell scheme. 

The expressions resulting from ( fO| ) can also be interpreted as loop graphs and are easily 
derived from the original Feynman graph. Applying again the FTT on the subtraction graphs, the 
resulting phase space integral over all tree graphs is UV finite. 

A closer inspection of the tree graphs of the FTT reveals, that the IR divergent parts solely arise 
from terms with a cut massless propagator. These can be immediately related to the corresponding 
IR divergent real emission graph by crossing the incoming piece of the cut propagator. 

Although in the limit of zero momentum-flow through this particle, the two contributions 
cancel each other, a suitable projection of the real emission part onto the virtual part is needed to 
cancel the IR peak locally in three-momentum space. For soft divergences, a simple approximation 
is obtained by setting the two contributions equal for momenta lower than a given soft energy E s . 

In the light of collinear divergences, a more sophisticated ansatz is to construct dipole sub- 
tractions, which cancel IR divergences locally also on the virtual side. The sum of the resulting 
analytic integrals for both parts of the dipoles would then be finite. This is left for future works. 



4. Threshold Singularities 



When the momentum integration is performed in the first line of ( |2.6| ), the integrand might get 
peaks in parts of the phase space where momenta of un-cut propagators are on-shell. The occur- 
rence of such peaks, although analytically integrable, leads to problems in the numerical evaluation 
of the integrand. In a series of papers for a direct numerical integration of massless amplitudes, 
Soper et al. ^ proposed a contour deformation of the integration into the complex plane to avoid 
the singularities. To allow for a direct implementation in an ordinary Monte Carlo generator, we 
add subtraction terms with zero real value but with the same peak structure to the integrand. These 
will cancel the peaks and allow for a better convergence in the numerical evaluation. 
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Full 1-Loop Corrections to QED Bhabha Scattering 
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Event Generation in Forwards Scattering Region, QED Bhabha Scattering 
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Figure 1: Differential cross sections for Monte Carlo integration and event generation of Bhabha scattering. 

In the sum of the tree graphs, peaks remain if pj t = (pi — pj) 2 > (m, + mj) 2 . The momentum 
constellation is such, that the two real particles could be produced at the same time. In the rest 
frame of p 2 h the peak of the threshold singularity is spherical with radius k s and a simple pole 



7(k)oc fk-k, 



%^{p% ,mf,m 2 j] 







(2|/#|) 



(4.1) 



X being the Kallen function. We can now evaluate the residue of the integrand and divide it by the 
singular structure. Adding the resulting function to the integrand symmetrically around the spher- 
ical peak, the peak will be cancelled and the integrated contribution of the fix function vanishes. 
More precisely, after a Lorentz transformation into the rest frame of p 2 t , the function we use reads 



4/4 



1 



kv 



2^ + ^J^ , )e(k'-(k,-c))0((k J +c, 



c- 



•k'). (4.2) 



This has to be transformed back to the integration frame and added to the integrand. Here, we also 
added a linear and cubic term in (k' — k s ), to make the joined integrand continuous and differen- 
tiable at the artificial borders introduced by the theta functions. 

It is possible, that for higher n-point functions there are overlapping threshold singularities. In 
this case, we can add further fix functions similar to ( |4~^ ) to cancel the peak structure. Since in this 
case, some of the added functions do not vanish after integration and give a small contribution to 



the final result, the width parameter c in (4.2) should be taken considerably small. This leads to a 
trade-off between accuracy and efficiency in the final integration. 



5. Results 

As a first application, we evaluate the one-loop cross section of Bhabha scattering in massive 
QED. The Bhabha scattering process is ideally suited to demonstrate the evaluation of processes 
at NLO by the Feynman Tree Theorem. The one-loop result for the cross section is well known 
and can easily be produced using automated loop-graph evaluation tools. Most complications 
inherent in the method - UV subtractions, IR cancellations, and threshold singularities - are present 
simultaneously. There are ten graphs in the one-loop corrections, which after rewriting them as tree 
graphs lead to a rich structure in the integrand. This has to be treated by multichannel integration 
methods. Furthermore, the smallness of the electron mass compared to the energy of 500 GeV 
where we evaluate the process provides a stringent test of the stability of the numerical integration. 
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We created analytical expressions for the loop graphs in computer-readable form using the 
Mathematica- and FORM-based packages FeynArts and FormCalc [0]. The matrix elements 
were then handed over to a specially crafted Mathematica program that creates subtraction graphs, 
cuts the loops and, where needed, calculates the fixing functions. For each tree graph, the program 
creates parameterizations (integration channels) which map the resulting phase space onto the unit 
hypercube, taking into account the peak structure. The resulting expressions for the matrix elements 
and the channels are then written out in Fortran code. 

As an integration routine, we choose the multi-channel algorithm VAMP [J8j|. We compare the 
final results for the cross section and angular distribution to an independent calculation that pro- 
ceeds along the usual way of integration via tensor reduction, using FeynArts and FormCalc. 
In the left plot of figure |l], we show the differential cross section for the full matrix element. The 
results are in complete agreement with FeynArts. The adaption of the grids in the multi-channel 
approach to the peaks works quite efficient. With 10000 sampling points, the error estimate on the 
numerical integration as returned by the Monte-Carlo integrator is less than 1%. 

The right hand plot of figure [j] shows results for event generation of the partonic QED process 
in the forwards scattering region. The distribution of events is again in agreement with predictions 
from FeynArts. Using the integrand obtained from the FTT, we generate unweighted events with 
a phase space including the on-shell momentum of a loop particle. This increases the dimension of 
integration from one to four. The efficiency of event generation is in the percent-regime. 

6. Conclusions 

We reported on a new method for computing NLO corrections to scattering cross sections. 
Here, all integrals are transformed into ordinary phase-space integrals (albeit with unusual bound- 
aries) that can be handled by an ordinary numerical multi-channel phase-space integrator. To this 
hand, we had not just to implement subtractions for UV and IR singularities, but furthermore sub- 
traction functions (fixing functions) for threshold singularities which do not cause problems in the 
usual semi-analytic methods. 

Extending the method to the full Standard Model and multi-particle processes, it promises 
important advantages over more conventional semi-analytic algorithms: The computational com- 
plexity does not increase dramatically with the number of legs in loop diagrams. Evaluating a NLO 
^-particle process should require similar CPU resources as a LO n+ 1 -particle process, summed 
inclusively over all particle species. Combining loop integration and phase-space sampling in a 
single step, we avoid a whole layer in the calculation. In particular, all terms are evaluated only up 
to the level of precision that is required by the actual simulation. 

There is still a long way before this method can actually improve the simulation of physics 
processes in the Standard Model or its extensions. On the one hand, we have to handle the more 
complicated IR behavior of QCD and state suitable renormalization conditions for the non-abelian 
theory. On the other hand, the method has to be augmented by a consistent treatment of unstable 
states, which appear in loops and, in our case, would become artificial external particles in event 
samples. Nevertheless, the method, has distinct advantages that warrant its further development 
towards realistic complete LHC and ILC applications. 
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